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We show in detail how the Jordan- Wigner transformation can be used to simulate any fermionic 
many-body Hamiltonian on a quantum computer. We develop an algorithm based on appropriate 
qubit gates that takes a general fermionic Hamiltonian, written as products of a given number of 
creation and annihilation operators, as input. To demonstrate the applicability of the algorithm, we 
calculate eigenvalues and eigenvectors of two model Hamiltonians, the well-known Hubbard model 
and a generalized pairing Hamiltonian. Extensions to other systems are discussed. 



I. INTRODUCTION 

A theoretical understanding of the behavior of many- 
body systems is a great challenge and provides fundamen- 
tal insights into quantum mechanical studies, as well as 
offering potential areas of applications. However, apart 
from some few analytically solvable problems, the typical 
absence of an exactly solvable contribution to the many- 
particle Hamiltonian means that we need reliable numer- 
ical many-body methods. These methods should allow 
for controlled approximations and provide a computa- 
tional scheme which accounts for successive many-body 
corrections in a systematic way. Typical examples of pop- 
ular many-body methods are coupled-cluster methods 
[U, 0, , various types of Monte Carlo methods [1, H, 0] , 
perturbative expansions 0, Q , Green's function methods 
the density-matrix renormalization group pd].[T3|. 
ab initio density functional theory Il3l . ll4l . llSlI and large- 
scale diagonalization methods [la, [3 ■ 

However, all these methods have to face in some form 
or the other the problem of an exponential growth in di- 
mensionality. For a system of P fermions which can be 
placed into N levels, the total number of basis states are 
f N\ 

given by I p I . The dimensional curse means that most 

quantum mechanical calculations on classical computers 
have exponential complexity and therefore are very hard 
to solve for larger systems. On the other hand, a so-called 
quantum computer, a particularly dedicated computer, 
can improve greatly on the size of systems that can be 
simulated, as foreseen by Feynman [S^, [2^. A quantum 
computer does not need an exponential amount of mem- 
ory to represent a quantum state. The basic unit of in- 
formation for a quantum computer is the so-called qubit 
or quantum bit. Any suitable two- level quantum system 
can be a qubit, but the standard model of quantum com- 
putation is a model where two-level quantum systems are 
located at different points in space, and are manipulated 
by a small universal set of operations. These operations 
are called gates in the same fashion as operations on bits 
in classical computers are called gates. 

For the example of P spin 1/2 particles, a classical 
computer needs 2^ bits to represent all possible states, 
while a quantum computer needs only P qubits. The 
complexity in number of qubits is thus linear. Based 



on these ideas, several groups have proposed various al- 
gorithms for simulating quantal many-body systems on 
quantum computers. Abrams and Lloyd, see for exam- 
ple Refs. [20l. [2H. introduced a quantum algorithm that 
uses the quantum fast Fourier transform to find eigen- 
values and eigenvectors of a given Hamiltonian, illustrat- 
ing how one could solve classically intractable problems 
with less than 100 qubits. Achieving a polynomial com- 
plexity in the number of operations needed to simulate a 
quantum system is not that straightforward however. To 
get efficient simulations in time one needs to transform 
the many-body Hamiltonian into a sum of operations on 
qubits, the building blocks of the quantum simulator and 
computer, so that the time evolution operator can be im- 
plemented in polynomial time. In Refs. 0,[2^[2^ it was 
shown how the Jordan- Wigner transformation in princi- 
ple does this for all Hamiltonians acting on fermionic 
many-body states. Based on this approach, recently two 
groups, see Refs. [27,, 23\^ published results where they 
used Nuclear Magnetic Resonance (NMR) qubits to sim- 
ulate the pairing Hamiltonian. 

The aim of this work is to develop an algorithm than 
allows one to perform a quantum computer simulation 
(or simply quantum simulation hereafter) of any many- 
body fermionic Hamiltonian. We show how to generate, 
via various Jordan- Wigner transformations, all qubit op- 
erations needed to simulate the time evolution operator 
of a given Hamiltonian. We also show that for a given 
term in an m-body fermionic Hamiltonian, the number of 
operations needed to simulate it is linear in the number 
of qubits or energy-levels of the system. The number of 
terms in the Hamiltonian is of the order of vn? for a gen- 
eral m-body interaction, making the simulation increas- 
ingly harder with higher order interactions. We special- 
ize our examples to a two-body Hamiltonian, since this 
is also the most general type of Hamiltonian encountered 
in many-body physics. Besides fields like nuclear physics, 
where three-body forces play a non-neglible role, a two- 
body Hamiltonian captures most of the relevant physics. 
The various transformations arc detailed in the next sec- 
tion. In Sec. Hh] we show in detail how to simulate a 
quantum computer finding the eigenvalues of any two- 
body Hamiltonian, with all available particle numbers, 
using the simulated time evolution operator. In that sec- 
tion we describe also the techniques which are necessary 
for the extraction of information using a phase-estimation 
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algorithm. 

To demonstrate the feasibihty of our algorithm, we 
present in Sec. |IV] selected results from applications of 
our algorithm to two simple modcl-Hamiltonians, a pair- 
ing Hamiltonian and the Hubbard model. We summarize 
our results and present future perspectives in Sec. |Vl 



II. ALGORITHM FOR QUANTUM 
COMPUTATIONS OF FERMIONIC SYSTEMS 

A. Hamiltonians 

A general two-body Hamiltonian for fermionic system 
can be written as 

H = Eo+YE^jalaj+ ^ V^jkiala^jUiak, (1) 

ij=l ijkl—1 

where Eq is a constant energy term, Eij represent all the 
one-particle terms, allowing for non-diagonal terms as 
well. The one-body term can represent a chosen single- 
particle potential, the kinetic energy or other more spe- 
cialized terms such as those discussed in connection with 
the Hubbard model [1^ or the pairing Hamiltonian dis- 
cussed below. The two-body interaction part is given by 
Vijki and can be any two-body interaction, from Coulomb 
interaction to the interaction between nucleons. The 
sums run over all possible single-particle levels N . Note 
that this model includes particle numbers from zero to 
the number of available quantum levels, n. To simulate 
states with fixed numbers of fermions one would have 
to either rewrite the Hamiltonian or generate specialized 
input states in the simulation. 

The algorithm which we will develop in this section and 
in Sec. IHII can treat any two-body Hamiltonian. How- 
ever, in our demonstrations of the quantum computing 
algorithm, we will limit ourselves to two simple models, 
which however capture much of the important physics in 
quantum mechanical many-body systems. We will also 
limit ourselves to spin j ~ 1/2 systems, although our al- 
gorithm can also simulate higher j-values, such as those 
which occur in nuclear, atomic and molecular physics, it 
simply uses one qubit for every available quantum state. 
These simple models are the Hubbard model and a pair- 
ing Hamiltonian. We start with the spin 1/2 Hubbard 
model, described by the following Hamiltonian 

i.a i,a 

+uJ2^U4-0■^-a^+, (2) 

i=l 

where and a are fermion creation and annihilation 
operators, respectively. This is a chain of sites where 
each site has room for one spin up fermion and one spin 
down fermion. The number of sites is iV, and the sums 
over a are sums over spin up and down only. Each site 



has a single-particle energy e. There is a repulsive term 
U if there is a pair of particles at the same site. It is 
energetically favourable to tunnel to neighbouring sites, 
described by the hopping terms with coupling constant 
-t. 

The second model-Hamiltonian is the simple pairing 
Hamiltonian 

Hp e^a\a, - ^5 X! al^ajaj , (3) 

The indices i and j run over the number of levels N, and 
the label z stands for a time-reversed state. The param- 
eter g is the strength of the pairing force while Si is the 
single-particle energy of level i. In our case we assume 
that the single-particle levels are equidistant (or degen- 
erate) with a fixed spacing d. Moreover, in our simple 
model, the degeneracy of the single-particle levels is set 
to 2j -I- 1 = 2, with j = 1/2 being the spin of the particle. 
This gives a set of single-particle states with the same 
spin projections as for the Hubbard model. Whereas in 
the Hubbard model we operate with different sites with 
spin up or spin down particles, our pairing models deals 
thus with levels with double degeneracy. Introducing the 
pair-creation operator = a|m^l-m' '^^^ rewrite 
the Hamiltonian in Eq. ^ as 

i ij>0 

where Ni = aja; is the number operator, and Si = id so 
that the single-particle orbitals are equally spaced at in- 
tervals d. The latter commutes with the Hamiltonian H. 
In this model, quantum numbers like seniority iS are good 
quantum numbers, and the eigenvalue problem can be 
rewritten in terms of blocks with good seniority. Loosely 
speaking, the seniority quantum number S is equal to 
the number of unpaired particles; see (30| for further de- 
tails. Furthermore, in a series of papers, Richardson, 
see for example Refs. [3l|, [H, H^, obtained the exact 
solution of the pairing Hamiltonian, with semi-analytic 
(since there is still the need for a numerical solution) ex- 
pressions for the eigenvalues and eigenvectors. The exact 
solutions have had important consequences for several 
fields, from Bose condensates to nuclear superconductiv- 
ity and is currently a very active field of studies, see for 
example Refs. [34l.[35j. Finally, for particle numbers up to 
P 20, the above model can be solved exactly through 
numerical diagonalization and one can obtain all eigen- 
values. It serves therefore also as an excellent ground 
for comparison with our algorithm based on models from 
quantum computing. 

B. Basic quantum gates 

Benioff showed that one could make a quantum me- 
chanical Turing machine by using various unitary opera- 
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tions on a quantum system, see Ref. [3^1 ■ Bcnioff demon- 
strated that a quantum computer can calculate anything 
a classical computer can. To do this one needs a quan- 
tum system and basic operations that can approximate 
all unitary operations on the chosen many-body system. 
We describe in this subsection the basic ingredients en- 
tering our algorithms. 

1. Qubits, gates and circuits 



Time^ 
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FIG. 1: A quantum circuit showing a single-qubit gate A and 
a two-qubit gate acting on a pair of qubits, represented by 
the horizontal lines. 



In this article we will use the standard model of quan- 
tum information, where the basic unit of information is 
the qubit, the quantum bit. As mentioned in the intro- 
duction, any suitable two-level quantum system can be 
a qubit, it is the smallest system there is with the least 
complex dynamics. Qubits are both abstract measures of 
information and physical objects. Actual physical qubits 
can be ions trapped in magnetic fields where lasers can 
access only two energy levels or the nuclear spins of some 
of the atoms in molecules accessed and manipulated by 
an NMR machine. Several other ideas have been pro- 
posed and some tested, see f37| . 

The computational basis for one qubit is |0) (repre- 
senting for example bit 0) for the first state and |1) (rep- 
resenting bit 1) for the second, and for a set of qubits the 
tensor products of these basis states for each qubit form 
a product basis. Below we write out the different basis 
states for a system of n qubits. 

|0) = I00---0) = |0) (g) |0) ® • • • (g) |0) 
|1) EE I00---1) = |0) ® |0) ® •••(g) |1) 



111' 



\1) \1) (E> ■ ■ ■ ® \1) 



(4) 



This is a 2"-dimensional system and we number the dif- 
ferent basis states using binary numbers corresponding 
to the order in which they appear in the tensor product. 

Quantum computing means to manipulate and mea- 
sure qubits in such a way that the results from a mea- 
surement yield the solutions to a given problem. The 
quantum operations we need to be able to perform our 
simulations are a small set of elementary singlc-qubit op- 
erations, or single-qubit gates, and one universal two- 
qubit gate, in our case the so-called CNOT gate defined 
below. 

To represent quantum computer algorithms graphi- 
cally we use circuit diagrams. In a circuit diagram each 
qubit is represented by a line, and operations on the dif- 
ferent qubits are represented by boxes. In fig. [1] we show 
an example of a quantum circuit, with the arrow indicat- 
ing the time evolution. The states |a) and \b) in the figure 
represent qubit states. In general, the total state will be 
a superposition of different qubit states. A single-qubit 
gate is an operation that only affects one physical qubit, 
for example one ion or one nuclear spin in a molecule. It 



is represented by a box on the line corresponding to the 
qubit in question. A single-qubit gate operates on one 
qubit and is therefore represented mathematically by a 
2x2 matrix while a two-qubit gate is represented by a 
4x4 matrix. As an example we can portray the so-called 
CNOT gate as matrix. 



(5) 



This is a very important gate, since one can show that 
it behaves as a universal two-qubit gate, and that we 
only need this two-qubit gate and a small set of single- 
qubit gates to be able to approximate any multi-qubit 
operation. One example of a universal set of single-qubit 
gates is given in Fig. [2l The products of these three 
operations on one qubit can approximate to an arbitrary 
precision any unitary operation on that qubit. 
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FIG. 2: Set of three elementary single-qubit gates and their 
matrix representations. The products of these three opera- 
tions on one qubit can approximate to an arbitrary precision 
any unitary operation on that qubit. 



2. Decomposing unitary operations into gates 

The next step is to find elementary operations on a set 
of qubits that can be put together in order to approxi- 
mate any unitary operation on the qubits. In this way 
we can perform computations on a quantum computer by 
performing many of these elementary operations in the 
correct order. 

There are three steps in finding the elementary oper- 
ations needed to simulate any unitary operation. First, 
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any dx d unitary matrix can be factorized into a product 
of at most d{d — l)/2 two-level unitary matrices, see for 
example Ref. [s^] for details. A two-level unitary ma- 
trix is a matrix that only acts non-trivially on two vector 
components when multiplied with a vector. For all other 
vector components it acts as the identity operation. 

The next step is to prove that any two-level unitary 
matrix can be implemented by one kind of two-qubit 
gate, for example the CNOT gate in Eq. (O, and single- 
qubit gates only. This simplifies the making of actual 
quantum computers as we only need one type of inter- 
action between pairs of qubits. All other operations are 
operations on one qubit at the time. 

Finally, these single-qubit operations can be approxi- 
mated to an arbitrary precision by a finite set of single- 
qubit gates. Such a set is called a universal set and one 
example is the phase gate, the so-called Hadamard gate 
and the 7r/8 gate. Fig. [2]shows these gates. By combining 
these properly with the CNOT gate one can approximate 
any unitary operation on a set of qubits. 



time evolution operator exp(— iTJAi) is found by using a 
Trotter approximation, for example 



U 



-iHAt 



j(5:,Hfc)At 



iH,.At 



(6) 

There are different ways to approximate U by products 
of exponentials of the different terms of the Hamiltonian, 
see Ref. [37[ and Eq. (|4ip . The essential idea is to find 
a form of the Hamiltonian where these factors in the ap- 
proximated time evolution operator can be further factor- 
ized into single- and two-qubit operations effectively. In 
Refs. [23, [11] it was shown how to do this in principle for 
any many-body fermion system using the Jordan- Wigner 
transformation. In this article we show how to create 
a quantum compiler that takes any many-body fermion 
Hamiltonian and outputs the quantum gates needed to 
simulate the time evolution operator. We implement 
it for the case of two-body fermion Hamiltonians and 
show results from numerical calculations finding the en- 
ergy levels of the well known pairing and Hubbard models. 



3. Quantum calculations 

The aspect of quantum computers we are focusing on 
in this article is their use in computing properties of 
quantum systems. When we want to use a quantum com- 
puter to find the energy levels of a quantum system or 
simulate it's dynamics, we need to simulate the time evo- 
lution operator of the Hamiltonian, U = exp{—iHAt). 
To do that on a quantum computer we must find a set 
of single- and two-qubit gates that would implement the 
time evolution on a set of qubits. For example, if we have 
one qubit in the state |a), we must find the single-qubit 
gates that when applied results in the qubit being in the 
state cxp{—iH At)\a) . 

From what we have written so far the naive way of 
simulating U would be to calculate it directly as a matrix 
in an appropriate basis, factorize it into two-level unitary 
matrices and then implement these by a set of universal 
gates. In a many-body fermion system for example, one 
could use the Fock basis to calculate C/ as a matrix. A 
fermion system with n different quantum levels can have 
from zero to n particles in each Fock basis state. A two- 
level system has four different basis states, |00), |01), |10) 
and 1 11), where |0) corresponds to an occupied quantum 
level. The time evolution matrix is then a 2" x 2" matrix. 
This matrix is then factorized into at most 2" (2" — l)/2 
two-level unitary matrices. An exponential amount of 
operations, in terms of the number of quantum levels, 
is needed to simulate U ; by definition not an effective 
simulation. 

This shows that quantum computers performing quan- 
tum simulations not necessarily fulfill their promise. For 
each physical system to be simulated one has to find rep- 
resentations of the Hamiltonian that leads to polynomial 
complexity in the number of operations. After one has 
found a proper representation of the Hamiltonian, the 



C. The Jordan- Wigner transformation 

For a spin-1/2 one-dimensional quantum spin-chain a 
fermionization procedure exists which allows the map- 
ping between spin operators and fermionic creation- 
annihilation operators. The algebra governing the spin 
chain is the SU{2) algebra, represented by the a- 
matrices. The Jordan- Wigner transformation is a trans- 
formation from fermionic annihilation and creation op- 
erators to the cr-matrices of a spin-1/2 chain, see for ex- 
ample Ref. [39I for more details on the Jordan- Wigner 
transformation. 

There is an isomorphism between the two systems, 
meaning that any a or operator can be transformed 
into a tensor product of cr-matrices operating on a set 
of qubits. This was explored by Somma et al. in 
Refs. [13, UHl- The authors demonstrated, with an em- 
phasis on single-particle fermionic operators, that the 
Jordan- Wigner transformation ensures efficient, i.e., not 
exponential complexity, simulations of a fermionic sys- 
tem on a quantum computer. Similar transformations 
must be found for other systems, in order to efficiently 
simulate many-body systems. This was the main point 
in Ref. [H. 

We present here the various ingredients needed in order 
to transform a given Hamiltonian into a practical form 
suitable for quantum mechanical simulations. 

We begin with the fermionic creation and annihilation 
operators, which satisfy the following anticommutation 
relations 

{ak,ai} ^ {al,aj} ^ 0, {al,ai}^5ki- (7) 

Thereafter we define the three traceless and Hermitian 
generators of the SU(2) group, the cr-matrices ax, a-y 
and cr^. Together with the identity matrix 1 they form 



5 



a complete basis for all Hcrmitian 2x2 matrices. They 
can be used to write all Hamiltonians on a spin 1/2 chain 
when taking sums of tensor products of these, in other 
words they form a product basis for the operators on the 
qubits. The three cr-matrices are 



1 

1 



-i 

1 



1 

-1 



(8) 



k*^ quantum level of the fermion system. When we have 
expressed the Hamiltonian as a sum of products of op- 
erations on the qubits representing the system, we must 
find a representation of the time evolution operator as 
products of two-qubit operations. These operations can 
be further decomposed into elementary operations, see 
subsection III B II for further discussion. 



We define the raising and lowering matrices as 



1 



1. 

2^ 



1 






1 



(9) 



The transformation is based on the fact that for each 
possible quantum state of the fermion system, there can 
be either one or zero fermions. Therefore we need n 
qubits for a system with n possible fermion states. A 
qubit in state |0)' = al\vacuum) represents a state with 
a fermion, while = \vacuum) represents no fermions. 
Then the raising operator (t_|- changes |1) into |0) when 



|0) 



(10) 



This means that a-^- acts as a creation operator, and 
tT_ acts as an annihilation operator. In addition, be- 
cause of the anticommutation of creation(annihilation) 
operators for different states we have ajajl vacuum) = 
—a2a\\vacuum) , meaning that for creation and annihi- 
lation operators for states higher than the state corre- 
sponding to the first qubit, we need to multiply with a 
CTz-matrix on all the qubits leading up to the one in ques- 
tion, in order to get the correct sign in the final opera- 
tion. This leads us to the Jordan- Wigner transformation 

MM 



n 



(11) 



The notation cr* cr;^ means a tensor product of the identity 
matrix on all qubits other than i and j, l(8)(Tz(8)lCg)cr+(g)l, 
Hi < j, with 1 being the identity matrices of appropriate 
dimension. 



D. Single-particle Hamiltonian 

What wc must do now is to apply the Jordan- Wigner 
transformation to a general fermionic Hamiltonian com- 
posed of creation and annihilation operators, so we can 
write it as a sum of products of a matrices. The matrix 
cr'^ is then an operation on the k^^ qubit representing the 



1. Jordan- Wigner transformation of the one-body part 

We first describe the procedure for the simplest case 
of a general single-particle Hamiltonian, 

i?i = ^ Enaja^ -\- ^ E,j (a^Oj + a]ai). (12) 

i i<j 

We now use the transformation of Eq. pTjl on the terms 

The diagonal terms of the one-particle Hamiltonian, 
that is the case where i = j, can be rewritten as 




2 ^ 



(13) 



since (cr^)^ — 1 which is the number operator. It counts 
whether or not a fermion is in state i. In the case of 
qubits counting whether or not the qubit is in state |0), 
we have eigenvalue one for |0) and eigenvalue zero for 
|1). The action of this Hamiltonian on qubit i can be 
simulated using the single-qubit operation 



U 



-i{l+cr^)Eii\t 



1 



(14) 



see subsection IIIB II for other examples of singlc-qubit 
gates. 

For the non-diagonal elements, i < j, not all of the 
az matrices multiply with each other and end up in the 
identity operation. As an example we will consider a five 
level system, n = 5, and look at the transformation of 
the term a|aj whith i = 2 and j = 4, 



0^04 = 1 (S 



SI cr+ 



51< 

1 Oz 



1 «) 1, 
D cr- 1, 



(15) 



The operation on all qubits before i and after j is identity, 
on qubits through j — 1 it is Gz ■ We can then write 
the non-diagonal one-body operators as 
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Using Eqs. (jTS]) and (fTB]) the total single-particle 
fermionic Hamiltonian of n quantum levels, transformed 
using the Jordan- Wigner transformation of Eq. (jlip . is 
written as 

i^i = ^ Eiia\ai + ^ Eij {a\aj + aja.;) 

i i<j 



i 




2. Transformation into two-qubit operations 

The Hamiltonian is now transformed into a sum of 
many-qubit operations, H = '^iHi- The a\ai term in 
Eq. (|15p for example, is transformed into a three-qubit 
operation. The next step is to factorize these many-qubit 
operations Hi into products of two-qubit operations, so 
that we in the end can get a product of two-qubit op- 
erations Ukij that when performed in order give us the 
time evolution operator corresponding to each term in 
the Hamiltonian, exp{—iHiAt) = Ylk ^ki- 

The first thing we do is to find a set of two-qubit op- 
erations that together give us the Hamiltonian, and later 
we will see that to find the time evolution from there is 
straightforward. The many-qubit terms in Eq. (|17|) are 
products of the type (TxCTz ■ ■ ■ o^Ux with or Oy at ei- 
ther end. These products have to be factorized into a 
series of two-qubit operations. What we wish to do is 
successively build up the operator using different unitary 
transformations. This can be achieved with successive 
operations with the a-matrices, starting with for example 
CT*, which can be transformed into tr!^, then transformed 
into a\^a'^^^ and so forth. Our goal now is to express 
each term in the Hamiltonian Eq. ([17]) as a product of 
the type < ct^ • • • ct^4 = (0^ Uw), with a dif- 

ferent form in the case where the Hamiltonian term starts 
and ends with a cr^ matrix. To achieve this we need the 



transformations in Eqs. (|A.8P - (jA.lip . We will use this to 
find the time-evolution operator for each Hamiltonian, 
see Eq. below. 

To understand how we factorize the Hamiltonian terms 
into single- and two-qubit operations we follow a bottom 
up procedure. First, if we have a two qubit system, with 
the operator az (8> 1, we see that the unitary operation 
exp(i7r/4(Tz ® cJz) transforms it into 



e-^'"'^''^®''^ {a, (g) 1) e'^/^"^®^^ =a,® a,. (18) 



In addition, if we start out with the operator cr* we 
can transform it into cr* or cr* using the operators 
exp(i7r/4cr,y) or exp(— Z7r/4cr^) accordingly. 

We can then generate the J^j, cr^' part of the terms 
in Eq. (jl7|) by succcsively applying the operator 
exp(i7r/4cr*cr^,) for / = 2 through I = j. Yielding the 
operator craY[k=i+i '^z with a phase of ±1, because of 
the sign change in Eqs. (jA.lOp and (jA.lip . We write ctq 
to show that we can start with both a ax and a ay ma- 
trix. To avoid the sign change we can simply use the 
operator exp{—iTT/4alal) instead for those cases where 
we have ay on site i instead of a^.. This way we always 
keep the same phase. 

Finally, we use the operator exp(i7r/4crj,) if we want 
the string of operators to end with ax, or exp(— i7r/4crj^) 
if we want it to end with ay. The string of opera- 
tors starts with either ax or ctj,. For an odd number 
of exp{zti7T / Aalal) operations, the operations that add 
a a^ to the string, the first operator has changed from 
what we started with. In other words we have ax in- 
stead of ay at the start of the string or vice versa, see 
Eqs. (jA.lOp and (jA.lip . By counting, we see that we do 
j — i of the exp(±i7r/4cr*cr^,) operations to get the string 
CT* cr*+^ • ■ ■ ai- and therefore if j — i is odd, the first matrix 
is the opposite of what we want in the final string. The 
following simple example can serve to clarify. We want 
the Hamiltonian a^a^a^ = ax®az®ax, and by using the 
transformations in Eqs. (|A.8P - (jA.ll|) we can construct it 
as a product of single- and two-qubit operations in the 
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following way, 



output state is the same as the input state. The operator 
term corresponding to Vijji has these equalities due to 
the anti-commutation relations 



^l<yl<yl- (19) 



We see that we have factorizcd crlalal into which means that 
Ul Ul Ul Ul a\ Ui U2 Us Ui . 

Now we can find the time-evolution operator 
exp(— iiJAt) corresponding to each term of the Hamilto- 
nian, which is the quantity of interest. Instead of start- 
ing with the operator tr* we start with the corresponding 
evolution operator and observe that 



a''jalajai 



(23) 



(24) 

The term aloja.iaj with i < j is described using the Pauli 
matrices 



(25) 



t/1'e-*'"--"C/ = (cos(a)l - i sin(a)cr2) U 
= cos(a)l — ism{a)U^azU 

_ -iU^cr^Ua 



(20) 



where a is a scalar. This means that we have a series 
of unitary transformations on this operator yielding the 
final evolution, namely 



(21) 



with the exact same unitary operations Uk as we find 
when we factorize the Hamiltonian. These are now the 
single- and two-qubit operations we were looking for, first 
we apply the operations Uk to the appropriate qubits, 
then exp(— icr^a) to qubit i, and then the ul operations, 
all in usual matrix multiplication order. 



E. Two-body Hamiltonian 

In this section we will do the same for the general two- 
body fcrmionic Hamiltonian. The two-body part of the 
Hamiltonian can be classified into diagonal elements and 
non-diagonal elements. Because of the Pauli principle 
and the anti-commutation relations for the creation and 
annihilation operators, some combinations of indices are 
not allowed. The two-body part of our Hamiltonian is 



= ^ Vi-jkiala^jUiak, 

ijkl 



(22) 



where the indices run over all possible states and n is the 
total number of available quantum states. The single- 
particle labels ijkl refer to their corresponding sets of 
quantum numbers, such as projection of total spin, num- 
ber of nodes in the single-particle wave function etc. 
Since every state ijkl is uniquely defined, we cannot have 
two equal creation or annihilation operators and there- 
fore i j and k ^l. 

When i = / and j = k, ov i = k and j = I, we have 
a diagonal element in the Hamiltonian matrix, and the 



\s=l J \t=l 
,t=l / \s=l 



16 



\t=i+l 



(26) 



When we add all four different permutations of i and j 
this is the number operator on qubit i multiplied with 
the number operator on qubit j. The eigenvalue is one if 
both qubits are in the state |0), that is the corresponding 
quantum states are both populated, and zero otherwise. 
We can in turn rewrite the sets of creation and annihila- 
tions in terms of the cr-matrices as 



(27) 



In the general case we can have three different sets 



of non-equal indices. Firstly, we see that ajaja/at; = 
al,a[ajai^ meaning that the exchange of i with k and j 



with / gives the same operator — ^ Vi 



ijkl 



Vkli 



This 



results in a two-body Hamiltonian with no equal indices 
HijM X! X! "^iokiiala^^aiak + a\a\ajai). (28) 

i<k j<l 

Choosing to order the indices from lowest to highest 
gives us the position where there will be (Tz-matrices to 
multiply with the different raising and lowering opera- 
tors, when we perform the Jordan- Wigner transforma- 
tion Eq. pip . The order of matrix multiplications is fixed 
once and for all, resulting in three different groups into 
which these terms fall, namely 



/ i < j <l < k, 

II i < I < j < k, 

III i < I < k < j, 



k I, 
k ^ I, 
k ^ I. 



(29) 
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These 12 possibilities for alojaiak are mirrored in the 
symmetric term in Eq. (j28|) giving us the 24 different 
possibihties when permuting four indices. 
The ijkl term of Eq. ((28]) is 

ajajazflfe + alajajUi = 

x(n'^0^-(n^^)^-- (30) 

In the case oii<j<l<kwe have 



n(a.)^) (n(-^)') 
n(-^)') (->v) (11' 



(31) 



Using Eq. (|A.12p . where we have the rules for sign 
changes when multiplying the raising and lowering op- 
erators with the cTz matrices, gives us 



„i „i „k-lk 



i i+1 1-1 j I l+l fe 

CT_cr^ ■ ■ ■ (ji (Tier I (T, ■■■a, a. 



(32) 



If we switch the order of i and j so that j < i < I < k, 
we change the order in which the cix-matrix is multiplied 
with the first raising and lowering matrices, resulting in 
a sign change. 

aja'jaiak + alajajai = 



.^i .^7 + 1 „l „k-lk 



(33) 



We get a change in sign for every permutation of the 
ordering of the indices from lowest to highest because 
of the matrix multiplication ordering. The ordering is 
described by another set of indices 

{Sa, 5/3,57,55} e {i,j,k,l} where < < < sg. 
We assign a number to each of the four indices, i 1, 



j <-!■ 2, / 3 and k A. If i < j < I < k we say 
the ordering isa = l,/3 = 2, 7 = 3 and (5 = 4, where 
a is a number from one to four indicating which of the 
indices i, j, I and k is the smallest. If i is the smallest, 
a = 1 and Sa = i- This allows us to give the sign of a 
given (ajajaittfc + aj^ajajai) term using the totally anti- 
symmetric tensor with four indices, which is -1-1 for even 
permutations, and —1 for odd permutations. For each 
of the three groups in Eq. (|29p we get a different set of 
raising and lowering operators on the lowest, next lowest 
and so on, indices, while the sign for the whole set is 
given by —£°'l^i^ . 

We are in the position where we can use the relation 
in Eq. ^ to express the Hamiltonian in terms of the a- 
matrices. We get 16 terms with products of four ax and 
or ay matrices in the first part of Eq. ()31|1 . then when we 
add the Hermitian conjugate we get another 16 terms. 
The terms with an odd number of ay matrices have an 
imaginary phase and are therefore cancelled out when 
adding the conjugates in Eq. ([28| . This leaves us with 
just the terms with four a^ matrices, four ay matrices 
and two of each in different orderings. The final result 
is given as an array with a global sign and factor given 
by the permutation of the ordering, and eight terms with 
different signs depending on which of the three groups, 
Eq. (|^^ . the set of indices belong to. These differing 
rules are due to the rules for cr^ multiplication with the 
raising and lowering operators, resulting in 



.t„t 



I II III 



+ 


+ 


+ 


0-5° o-z" 


azax 


a'^x Oz 


■■■Oz 


0%' 




+ 


+ 




■ ■ Ox 


ay ■ ■ 


ay 




+ 




+ 




■■ay 


Ox ■ ■ 


ay 




< + 


+ 




<7x 


■■ay 


ay ■ ■ 


Ox 


(34) 


+ 


+ 




ay 


■■Ox 


Ox ■ ■ 


ay 




+ 




+ 


ay 


■ ■ Ox 


ay ■ ■ 


Ox 






+ 


+ 


ay 


■■ay 


Ox ■ ■ 


Ox 






+ 


+ 


ay 


■■ay 


ay ■ ■ 


ay 





where the letters /, // and /// refer to the subgroups 

defined in Eq. (gH). 

As for the single-particle operators in subsection III Dl 
we now need to factorize these multi-qubit terms in the 
Hamiltonian to products of two-qubit and single-qubit 
operators. Instead of transforming a product of the form 
az - ■ ■ zb, we now need to transform a product of the form 
az ■ ■ ■ zbcz ■ ■ ■ zd, where a, &, c and d are short for either 
ax or ay while z is short for CTz. The generalization is 
quite straightforward, as we see that if the initial operator 
is (t|°(7^"' instead of just erf" , we can use the same set of 
transformations as for the single-particle case, 

c/^••c/lV^c/l•••c/fc 



= CTq CTz ■ • • OzOf, 

Ul^^uUl-al^Ui^ 
= ol^Oz ■ ■■OzO^^'al 



■Uk 



(35) 
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Using the same unitary two-qubit transformations, which 
we now call V, that take crT' to (Tc^cTz • • • cTz^'d' ^ ^'^'^ 

Vs--- V^Ul ■ ■ ■ Ulal-al- Ui---UkVi--- V, 

= cr^°fT^ • • • Tzcr^cr^-'cr, ■ ■ ■ a^a^ (36) 

This straightforward generalization of the procedure from 
the single-particle Hamiltonian case is possible because 
the operations commute when performed on different 
qubits. 

With the above expressions, we can start with the uni- 
tary operator exp(— zacrf^crl"') and have two different se- 
ries of unitary operators that give us the evolution op- 
erator of the desired Hamiltonian. The U operators are 
defined as in Eq. (PT|) . 

while the V operators are defined in a similar way 

where the (T-matrices without subscripts represent that 
we can have <Jx or Uy in each position. 

This gives us the total evolution operator for each term 
in Eq. ^ 




Here we have all the single- and two-qubit operations 
we need to perform on our set of qubits, that were 
initially in the state \ip), to simulate the time evolu- 
tion exp(—i_fffe At) !■(/;) of the Hamiltonian term Hk = 
a'^^az ■ ■ ■ OzO^f' a^'i Oz ■ ■ ■ <Jz<^^^ ■ Every factor in the above 
equation is a single- or two-qubit operation that must be 
performed on the qubits in proper matrix multiplication 
order. 

When using the Jordan- Wigner transformation of 
Eq. (|lip applied to our two model Hamiltonians of 
Eqs. (HI and ([3]), we choose a representation with two 
qubits at each site. These correspond to fermions with 
spin up and down, respectively. The number of qubits, 
n. is always the total number of available quantum states 
and therefore it is straightforward to use this model on 
systems with higher degeneracy, such as those encoun- 
tered in quantum chemistry or nuclear physics 
Site one spin up is qubit one, site one spin down is qubit 
two and site two spin up is qubit three and so on. To 
get all the quantum gates one needs to simulate a given 
Hamiltonian one needs to input the correct Eij and Vijki 
values. 



F. Complexity of the quantum computing 
algorithm 

In order to test the efficiency of a quantum algorithm, 
one needs to know how many qubits, and how many oper- 
ations on these, are needed to implement the algorithm. 
Usually this is a function of the dimension of the Hilbert 
space on which the Hamiltonian acts. The natural in- 
put scale in the fermionic simulator is the number of 
quantum states, n, that are available to the fermions. 
In our simulations of the Hubbard and the pairing mod- 
els of Eqs. (HI and ([3]), respectively, the number of qubits 
is n = IN since we have chosen systems with double- 
degeneracy for every single-particle state, where N is 
the number of energy- levels in the model. We use one 
qubit to represent each possible fermion state, on a real 
quantum computer, however, one should implement some 
error-correction procedure using several qubits for each 
state, see Ref. [37| . The complexity in number of qubits 
remains linear, however, since Oin) qubits are needed for 
error correction. 

The single- particle Hamiltonian has potentially 0(r?^ 
different terms. The two-particle Hamiltonian has up 
to 0(n^^ Vijki terms. A general m-body interaction has 
in the worst case 0{n^"^) terms. It is straightforward to 
convince oneself that the pairing model has 0{n^) terms 
while in the Hubbard model we end up with 0(n) terms. 
Not all models have maximum complexity in the different 
m-body interactions. 

How many two-qubit operations do each of these terms 
need to be simulated? First of all a two-qubit opera- 
tion will in general have to be decomposed into a series 
of universal single- and two-qubit operations, depending 
entirely on the given quantum simulator. A particular 
physical realization might have a natural implementa- 
tion of the al (g) ai gate and save a lot of intermediary 
operations. Others will have to use a fixed number of 
operations in order to apply the operation on any two 
qubits. A system with only nearest neighbor interactions 
would have to use 0{n) operations for each a* (8) cr^ gate, 
and thereby increase the polynomial complexity by one 
degree. 

In our discussion on the one-body part of the Hamilto- 
nian, we saw that for each Eij we obtained the a|aj -l-a|ai 
operator which is transformed into the two terms in 
Eq. (jl6p . (7x(Jz ■ ■ ■ o'zO'x and UyCTz ■ ■ ■ Uz<^y Wc showed 
how these terms are decomposed into j~i + 2 operations, 
leading to twice as many unitary transformations on an 
operator, VAV^ for the time evolution. The average of 
j — i is rt/2 in this case and in total we need to perform 
2 X 2 X n/2 = 2n two-qubit operations per single-particle 
term in the Hamiltonian, a linear complexity. 

In the two-particle case each term T^yfci(a|a|a/afc -I- 
aj^ajajai) is transformed into a sum of eight operators 
of the form a-'^az ■ ■ ■ cTzfr'^f a"-' az ■ ■ ■ cr^cr■'^ Eq. The 
two parts of these operators are implemented in the same 
way as the a^az ■ ■ ■ o'zO'-' term of the single-particle Hamil- 
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tonian, which means they require — and sg — s-y op- 
erations, since Sa < S/3 < < ss the average is n/4. For 
both of these parts we need to perform both the unitary 
operation V and it's Hermitian conjugate V'' . In the end 
we need 2x2x8xn/4 = 8n two-qubit operations per 
two-particle term in the Hamiltonian, the complexity is 
linear. 

A term of an m-body Hamiltonian will be transformed 
into 2^™ operators since each annihilation and creation 
operator is transformed into a sum of ax and <7y matri- 
ces. All the imaginary terms cancel out and we are left 
with 2^™"^ terms. Each of these terms will include 2m 
a matrices, in products of the form OfcLi (^^(^z ' ' ' ^zcr^ , 
and we use the same procedure as discussed above to 
decompose these m factors into unitary transformations. 
In this case each factor will require an average of n/2m 
operations for the same reasons as in the two-body case. 
All in all, each m-body term in the Hamiltonian requires 
22m-i X 2 X m X n/2m = 2^'"~^n operations. 

Thus, the complexity for simulating one m-body term 
of a fermionic many-body Hamiltonian is linear in the 
number of two-qubit operations, but the number of terms 
is not. For a full-fledged simulation of general three- 
body forces, in common use in nuclear physics [iol . l4ll . 
l42j . the total complexity of the simulation is 0{'n7). A 
complete two-particle Hamiltonian would be C'(n^).The 
bottleneck in these simulations is the number of terms 
in the Hamiltonian, and for systems with less than the 
full number of terms the simulation will be faster. This 
is much better than the exponential complexity of most 
simulations on classical computers. 



III. ALGORITHMIC DETAILS 

Having detailed how a general Hamiltonian, of two- 
body nature in our case, can be decomposed in terms of 
various quantum gates, we present here details of the im- 
plementation of our algorithm for finding eigenvalues and 
eigenvectors of a many-fermion system. For our tests of 
the fermionic simulation algorithm we have implemented 
the phase-estimation algorithm from [s^l which finds the 
eigenvalues of an Hamiltonian operating on a set of sim- 
ulation qubits. There are also other quantum computer 
algorithms for finding expectation values and correlation 
functions, as discussed by Somma et al. in Refs. [25l.l26|. 
In the following we first describe the phase-estimation al- 
gorithm, and then describe its implementation and meth- 
ods we have developed in using this algorithm. A much 
more thorough description of quantum computers and 
the phase-estimation algorithm can be found in [4^. 

A. Phase-estimation algorithm 

To find the eigenvalues of the Hamiltonian we use the 
unitary time evolution operator we get from the Hamil- 
tonian. We have a set of simulation qubits representing 



the system governed by the Hamiltonian, and a set of 
auxiliary qubits, called work qubits [13, [2l|, in which we 
will store the eigenvalues of the time evolution operator. 
The procedure is to perform several controlled time evo- 
lutions with work qubits as control qubits and the simu- 
lation qubits as targets, see for example Ref. [13] for in- 
formation on controlled qubit operations. For each work 
qubit we perform the controlled operation on the simu- 
lation qubits with a different time parameter, giving all 
the work qubits different phases. The information stored 
in their phases is extracted using first an inverse Fourier 
transform on the work qubits alone, and then performing 
a measurement on them. The values of the measurements 
give us directly the eigenvalues of the Hamiltonian after 
the algorithm has been performed a number of times. 

The input state of the simulation qubits is a ran- 
dom state in our implementation, which is also a ran- 
dom superposition of the eigenvectors of the Hamilto- 
nian = y^t. C k\k). It does not have to be a random 
state, and iri[4J| the authors describe a quasi-adiabatic 
approach, where the initial state is created by starting 
in the ground state for the non-interacting Hamiltonian, 
a qubit basis state, e.g. jOlOl • • • 101), and then slowly 
the interacting part of the Hamiltonian is turned on. 
This gives us an initial state mostly comprising the true 
ground state, but it can also have parts of the lower ex- 
cited states if the interacting Hamiltonian is turned on 
a bit faster. In for example nuclear physics it is com- 
mon to use a starting state for large-scale diagonaliza- 
tions that reflects some of the features of the states one 
wishes to study. A typical example is to include pairing 
correlations in the trial wave function, see for example 
Refs. [3 US- Iterative methods such as the Lanczo's di- 
agonalization technique [l^ . l45j converge much faster if 
such starting vectors are used. However, although more 
iterations are needed, even a random starting vector con- 
verges to the wanted states. 

The final state of all the qubits after an inverse Fourier 
transform on the work qubits is 

^Cfc|0W2*)®|A:). (40) 

k 

If the algorithm works perfectly, \k) should be an exact 
eigenstate of U , with an exact eigenvalue 4>^^^ ■ When we 
have the eigenvalues of the time evolution operator we 
easily find the eigenvalues of the Hamiltonian. We can 
summarize schematically the phase-estimation algorithm 
as follows: 

1. Intialize each of the work qubits to 1/V2(|0) -f- |1)) 
by initializing to |0) and applying the Hadamard 
gate, H, see Fig. [2] 

2. Initialize the simulation qubits to a random or spec- 
ified state, depending on the whether one wants the 
whole eigenvalue spectrum. 

3. Perform conditional time evolutions on the simula- 
tion qubits, with different timesteps At and differ- 
ent work qubits as the control qubits. 
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|0) +e2"'(2™ '-^'ll) = |0) +e2"'''-'^-|: 
jO) +e2"'(2''^'il) = \Q) ^ e^^^o.4,2-<i>^. 
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FIG. 3: Phase estimation circuit showing all the different qubit lines schematically with operations represented by boxes. The 
boxes connected by vertical lines to other qubit lines are controlled operations, with the qubit with the black dot as the control 
qubit. 



4. Perform an inverse Fourier transform on the work 
qubits. 

5. Measure the work qubits to extract the phase. 

6. Repeat steps 1-6 until the probabihty distribution 
I 



gathered from the measurement results is good 
enough to read out the wanted eigenvalues. 



As discussed above a set of two-qubit operations can 
be simulated by the CNOT two-qubit operation and a 
universal set of singlc-qubit operations. We will not use 
or discuss any such implementation in this article, as one 
will have to use a different set for each physical realiza- 
tion of a quantum computer. When simulating a fermion 
system with a given quantum computer, our algorithm 
will first take the fermionic many-body evolution opera- 
tor to a series of two-qubit and single-qubit operations, 
and then one will have to have a system dependent setup 
that takes these operations to the basic building blocks 
that form the appropriate universal set. 



In subsection III El we showed how to take any two- 
particle fermionic Hamiltonian to a set of two-qubit op- 
erations that approximate the evolution operator. In ad- 
dition we must use one of the Trotter approximations 
[H, 113, Sil Eqs. dH) and (gl]) that take the evolution op- 
erator of a sum of terms to the product of the evolution 
operator of the individual terms, see for example Ref. [s^l 
for details. To order C(At^) in the error we use Eq. ^ 
while to order 0{/S.t^) we have 



-i{A+B)At 



-iAAt/2 -iBAt -iAAt/2 



0{At^). (41) 



B. Output of the phase-estimation algorithm 

The output of the phase-estimation algorithm is a se- 
ries of measurements of the w number of work qubits. 
Putting them all together we get a probability distribu- 
tion that estimates the amplitudes |cfep for each eigen- 
value The 0['=l2"' values we measure from the work 
qubits, see Eq. (|40p . are binary numbers from zero to 
2"' — 1, where each one translates to a given eigenvalue of 
the Hamiltonian depending on the parameters we have 
used in our simulation. When accurate, a set of simu- 
lated measurements will give a distribution with peaks 
around the true eigenvalues. The probability distribu- 
tion is calculated by applying non-normalized projection 
operators to the qubit state, 

(|^W2*)((/)W2*|®l) |^^c,|0,2*)^|i)j =Cfc|</)W2*)®|fc). 

The length of this vector squared gives us the probability, 

Cfc|<^W)2* ® |fc)|' = |c,p(0W2*|(/)W2*)(fc|fc) = Icfep. 

(42) 

Since we do not employ the exact evolution due to differ- 
ent approximations, we can have non-zero probabilities 
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for all values of yielding a distribution without sharp 
peaks for the correct eigenvalues and possibly peaks in 
the wrong places. If we use different random input states 
for every run through the quantum computer and gather 
all the measurements in one probability distribution, all 
the eigenvectors in the input state = Ck\k) should 
average out to the same amplitude. This means that 
eigenvalues with higher multiplicity, i.e., higher degen- 
eracy, will show up as taller peaks in the probability 
distribution, while non-degenerate eigenvalues might be 
difficult to find. 

To properly estimate the eigenvalues Ek of the Hamil- 
tonian from this distribution, one must take into account 
the periodicity of e^'^''^. If < 0' < 1 and = 0' + s, 
where s is an integer, then e^'^*'^ = e^'^"^ . This means 
that to get all the eigenvalues correctly (j) must be positive 
and less than one. Since (j) = —EkA-t/2iT this means all 
the eigenvalues Ek must be negative, this merely means 
subtracting a constant we denote E„^ax from the Hamil- 
tonian, H' = H — Emax, where Emax is greater than the 
largest eigenvalue of H. The values we read out from 
the work qubits are integers from zero to 2™ — 1. In 
other words, we have 0['=l2"' e [0, 2^ - 1], with = for 
At = 0. 

The value = corresponds to the lowest eigen- 
value possible to measure, Emin, while = 1 corre- 
sponds to Emax- The interval of possible values is then 
Emax — Emin = 27r/Ai. If wc waut to have all possible 
eigenvalues in the interval the largest value Ai can have 
is 

27r 

max(At) = — (43) 



1. Spectrum analysis 

In the general case one does not know the upper and 
lower bounds on the eigenvalues beforehand, and there- 
fore for a given Emax and At one does not know if the 
0[''l values are the correct ones, or if an integer has been 
lost in the exponential function. 

When = 0' + s for one At, and we slightly change A^, 
0' will change if s as the period of the exponential 
function is a function of At. To find out which of 0['^ls are 
greater than one, we perform the phase-estimation algo- 
rithm with different values for At and see which eigen- 
values shift. If we measure the same after adding 6t to 
the time step, and {At + 5t)/At is not a rational number, 
we know that < 1. In practice it docs not have to be 
an irrational number, but only some unlikely fraction. 

There are at least two methods for finding the eigen- 
values. One can start with a large positive Emax and a 
small At, hoping to find that the whole spectrum falls 
within the range [Emin, Emax]i and from there zoom in 
until the maximal eigenvalue is slightly less than Emax 
and the groundstate energy is slightly larger than Emm- 
This way the whole spectrum is covered at once. From 



there we can also zoom in on specific areas of the spec- 
trum, searching the location of the true eigenvalues by 
shifting At. 

The number of measurements needed will depend en- 
tirely on the statistics of the probability distribution. 
The number of eigenvalues within the given energy range 
determines the resolution needed. That said, the number 
of measurements is not a bottleneck in quantum com- 
puter calculations. The quantum computer will prepare 
the states, apply all the operations in the circuit and 
measure. Then it will do it all again. Each measure- 
ment will be independent of the others as the system is 
restarted each time. This way the serious problem of de- 
coherence only apply within each run, and the number 
of measurements is only limited by the patience of the 
scientists operating the quantum computer. 



IV. RESULTS AND DISCUSSION 

In this section we present the results for the Hubbard 
model and the pairing model of Eqs. ^ and 13]), respec- 
tively, and compare the simulations to exact diagonaliza- 
tion results. In Fig. |4] we see the resulting probability 
distribution from the simulated measurements, giving us 
the eigenvalues of the pairing model with six degenerate 
energy levels and from zero to 12 particles. The pairing 
strength was set to g = 1. The eigenvalues from the exact 
solutions of these many-particle problems are 0, —1, —2, 
—3, —4, —5, —6, —8, —9, —12. All the eigenvalues are 
not seen as this is the probability distribution resulting 
from one random input state. A different random in- 
put state in each run could be implemented on an actual 
quantum computer. These are results for the degenerate 
model, where the single-particle energies of the doubly 
degenerate levels are set to zero for illustrate purposes 
only, since analytic formula are available for the exact 
eigenvalues. The block diagonal structure of the pairing 
Hamiltonian has not been used to our advantage in this 
straightforward simulation as the qubit basis includes all 
particle numbers. 

We have also performed tests of the algorithm for the 
non-degenerate case, with excellent agreement with our 
diagonalization codes, see discussion in Ref. [131. This is 
seen in Fig. [5] where we have simulated the pairing model 
with four energy levels with a total possibility of eight 
fermions. We have chosen .g = 1 and d = 0.5, so this is a 
model with low degcneray and since the dimension of the 
system is 2^ = 256 there is a lot of different eigenvalues. 
To find the whole spectrum one would have to employ 
some of the techniques discussed in subsection IIIIBI 



A. Number of work qubits versus number of 
simulation qubits 

The largest possible amount of different eigenvalues is 
2*^, where s is the number of simulation qubits. The 
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FIG. 4: Resulting probability distribution from the simulated 
measurements, giving us the eigenvalues of the pairing model 
with six degenerate energy levels with a total possibility of 
12 particles and pairing strength g = 1. The correct eigen- 
values are 0, -1, -2, -3, -4, -5, -6, -8, -9, -12. All the 
eigenvalues are not seen as this is the probability distribution 
resulting from one random input state. A different random 
input state in each run could be implemented on an actual 
quantum computer and would eventually yield peaks of height 
corresponding to the degeneracy of each eigenvalue. 
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chose t = was just because of the higher degeneracy and 
therefore fewer eigenvalues. The number of work qubits 
is 16 and the number of simulation qubits is eight for a 
total of 24 qubits. The difference between work qubits 
and simulation qubits is eight which means there are 2* 
possible energy values for each eigenvalue. Combining 
that with the high degeneracy we get a very sharp reso- 
lution. The correct eigenvalues with degeneracies arc ob- 
tained from exact diagonalization of the Hamiltonian, the 
degeneracy follows the eigenvalue in paranthesis: 0(1), 
1(8), 2(24), 3(36), 4(40), 5(48), 6(38), 7(24), 8(24), 9(4), 
10(8), 12(1). We can clearly see that even though we 
have a random input state, with a random superposition 
of the eigenvectors, there is a correspondence between 
the height of the peaks in the plot and the degeneracy of 
the eigenvalues they represent. 




FIG. 6: The energy levels of the Hubbard model of Eq. ((2|, 
simulated with a total of 24 qubits, of which eight were sim- 
ulation qubits and 16 were work qubits. In this run we chose 
e = 1, t = and U = 1. The reason we chose t = was just 
because of the higher degeneracy and therefore fewer eigenval- 
ues. The correct eigenvalues are obtained from exact diago- 
nalization, with the level of degeneracy following in paranthe- 
sis: 0(1), 1(8), 2(24), 3(36), 4(40), 5(48), 6(38), 7(24), 8(24), 
9(4), 10(8), 12(1). 



FIG. 5: The eigenvalues of the non-degenerate pairing model 
with four energy levels with a total possibility of 8 particles, 
the level spacing d is 0.5 and the pairing strength g is 1. The 
correct eigenvalues are obtained from exact diagonalization, 
but in this case there is a multitude of eigenvalues and only 
some eigenvalues are found from this first simulation. 



resolution in the energy spectrum we get from measuring 
upon the work qubits is 2"", with w the number of work 
qubits. Therefore the resolution per eigenvalue in a non- 
degenerate system is 2™"*. The higher the degeneracy 
the less work qubits are needed. 

In Fig. [H] we see the results for the Hubbard model 
Eq. ^ with e = 1, t = and U — 1. The reason we 



B. Number of time intervals 

The number of time intervals, /, is the number of 
times we must implement the time evolution operator 
in order to reduce the error in the Trotter approximation 
[H, [13, [13 , ss*^ ^'i- ® • "^^1' program wc have only 
implemented the simplest Trotter approximation and in 
our case wc find that we do not need a large / before 
the error is small enough. In Fig. [S] / is only one, but 
here we have a large number of work qubits. For other 
or larger systems it might pay off to use a higher order 
Trotter approximation. The total number of operations 
that have to be done is a multiple of /, but this number 



14 



also increases for higher order Trotter approximations, so 
for each case there is an optimal choice of approximation. 

In Figs. [7] and [5] we see the errors deriving from the 
Trotter approximation, and how they are reduced by in- 
creasing the number of time intervals. The results in this 
figure are for the degenerate pairing model with 24 qubits 
in total, and ten simulation qubits with d = and g ~ 1- 
In Fig. [7] we had 1=1 while in Fig. [5] / was set to ten. 
Both simulations used the same starting state. The er- 
rors are seen as the small spikes around the large ones 
which represent some of the eigenvalues of the system. 
The exact eigenvalues are 0, —1, —2, —3, —4, —5, —6, 
-8, -9. 
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FIG. 7: Pairing model simulated with 24 qubits, where 14 
were simulation qubits, i.e. there are 14 available quantum 
levels, and 10 were work qubits. The correct eigenvalues are 
0, —1, —2, —3, —4, —5, —6, —8, —9. In this run we did not 
divide up the time interval to reduce the error in the Trotter 
approximation, i.e., 7 = 1. 



C. Number of operations 

Counting the number of single-qubit and cr^cr^ opera- 
tions for different sizes of systems simulated gives us an 
indication of the decoherence time needed for different 
physical realizations of a quantum simulator or computer. 
The decoherence time is an average time in which the 
state of the qubits will be destroyed by noise, also called 
decoherence, while the operation time is the average time 
an operation takes to perform on the given system. Their 
fraction is the number of operations possible to perform 
before decoherence destroys the computation. In table |I] 
we have listed the number of gates used for the pairing 
model, Hp, and the Hubbard model, Hh, for different 
number of simulation qubits. 
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FIG. 8: Pairing model simulated with 24 qubits, where 14 
were simulation qubits, i.e. there are 14 available quantum 
levels, and 10 were work qubits. The correct eigenvalues are 
0, —1, —2, —3, —4, —5, —6, —8, —9. In this run we divided 
the time interval into 10 equally space parts in order to reduce 
the error in the Trotter approximation, i.e., / = 10. 





s = 2 s = 4 s = 6 s = 8 


s = 10 s = 12 


Hp 


9 119 333 651 


1073 1598 


Hh 


9 51 93 135 


177 219 



TABLE I: Number of two-qubit gates used in simulating 
the time evolution operator of the pairing model. Hp, and 
the Hubbard model, Hh, for different number of simulation 
qubits s. 



V. CONCLUSION 

In this article we have shown explicitly how the Jordan- 
Wigner transformation is used to simulate any many- 
body fermionic Hamiltonian by two-qubit operations. 
We have shown how the simulation of such Hamiltonian 
terms of products of creation and annihilation operators 
are represented by a number of operations linear in the 
number of qubits. To perform efficient quantum simula- 
tions on quantum computers one needs transformations 
that take the Hamiltonian in question to a set of opera- 
tions on the qubits simulating the physical system. An 
example of such a transformation employed in ths work, 
is the Jordan- Wigner transformation. With the appro- 
priate transformation and relevant gates or quantum cir- 
cuits, one can taylor an actual quantum computer to sim- 
ulate and solve the eigenvalue and eigenvector problems 
for different quantum systems. Specialized quantum sim- 
ulators might be more efficient in solving some problems 
than others because of similarities in algebras between 
physical system of qubits and the physical system simu- 
lated. 

We have limited the applications to two simple and 
well-studied models that provide, via exact eigenvalues. 
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a good testing ground for our quantum computing based 
algorithm. For both the pairing model and the Hubbard 
model we obtain an excellent agreement. We plan to 
extend the area of application to quantum mechanical 
studies of systems in nuclear physics, such as a compari- 
son of shell-model studies of oxygen or calcium isotopes 
where the nucleons are active in a given number of single- 
particle orbits 0, [ill • These single-particle orbits have 
normally a higher degeneracy than 2, a degeneracy stud- 
ied here. However, the algorithm we have developed al- 
lows for the inclusion of any degeneracy, meaning in turn 
that with a given interaction Vijki and single-particle en- 
ergies, we can compare the nuclear shell-model (configu- 
ration interaction) calculations with our algorithm. 
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APPENDIX: USEFUL RELATIONS 

We list here some useful relations involving different a 
matrices, 

CTxCTz = -io-y, azCTx = icTy, [ax,crz] = -2iay, (A.l) 

(JxCTy = i(Jz, Cry(7x = -iCTz, ['Jx,(7y] = 2i(Tz, (A. 2) 

and 

ffyffz = icTx, CFzCTy = -icTx, Wy^CTz] = 2i(Tj;. (A. 3) 

For any two non-equal cr-matrices a and b we have 

aba = -b. (A. 4) 



The Hcrmitian cr-matrices cTx, cry and result in the 
identity matrix when squared 

^2 = 1, al = l, al = X (A.5) 

which can be used to obtain simplified expressions for 
exponential functions involving cr-matrices 

^±^aa ^ cos(a)l ± i sin(a)cr. (A.6) 

The equations we list below are necessary for the re- 
lation between a general unitary transformation on a 
set of qubits with a product of two-qubit unitary trans- 
formations. We have the general equation for a, 5 G 
{(Jx.Uy, cr^}, where a^b. 

e-'"/4'^6e"/4« = i (1 - ia)h[\ + la) 

= - (5 -f a6a -I- z [6, a] ) 

= '-[b,a]. (A.7) 
The more specialized equations read 



e~ 


'i-K/ia^ 'iir/4cr^. _ 


(A.8) 


e~ 




(A.9) 


e' 


'^^/'^^^axe'"/^"' =cjy, 


(A.IO) 


e' 




(A.ll) 



We need also different products of the operatorcr^ with 
the raising and lowering operators 



(J^Oz = -cr+ 


(A.12) 


OzCTj^ = cr+. 


(A.13) 


a^az = cr_. 


(A. 14) 


(7zO'- = — cr_. 


(A.15) 




(A.16) 
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